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Being a true two-dimensional crystal, graphene has special properties. In particular, a point-like 
defect in graphene may induce perturbations in the long range. This characteristic questions the 
validity of using a supercell geometry in an attempt to explore the properties of an isolated defect. 
Still, this approach is often used in ab-initio electronic structure calculations, for instance. How 
does this approach converge with the size of the supercell is generally not tackled for the obvious 
reason of keeping the computational load to an affordable level. The present paper addresses the 
problem of substitutional nitrogen doping of graphene. DFT calculations have been performed 
for 9x9 and 10 x 10 supercells. Although these calculations correspond to N concentrations 
that differ by ~ 10%, the local densities of states on and around the defects are found to depend 
significantly on the supercell size. Fitting the DFT results by a tight-binding Hamiltonian makes 
it possible to explore the effects of a random distribution of the substitutional N atoms, in the case 
of finite concentrations, and to approach the case of an isolated impurity when the concentration 
vanishes. The tight-binding Hamiltonian is used to calculate the STM image of graphene around an 
isolated N atom. STM images are also calculated for graphene doped with 0.5 at% concentration of 
nitrogen. The results are discussed in the light of recent experimental data and the conclusions of 
the calculations are extended to other point defects in graphene. 

PACS numbers: 73.22.Pr, 31.15.aq 
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I. INTRODUCTION 

Local defects and chemical doping is a well- 
documented way to tune the electronic properties of 
graphenei More specifically, the benefits of doping have 
been underlined in the context of (bio)sensing; 2 lithium 
incorporation battery^^ and in other fields. Nitrogen 
(N) is a natural substitute for carbon in the honeycomb 
structure due to both its ability to form sp 2 bonds and 
its pentavalent character. It is not a surprise, then, that 
many publications deal with the production of N-doped 
graphene realized either by direct growth of modified lay- 
ers^ or by post-synthesis treatments. 7,8 Doping single- 
wall carbon nanotubes has also been considered^— 

Recent STM-STS experiments of N-doped graphene^ 
have demonstrated the occurrence of chemical substitu- 
tion. These experiments provide us with a detailed and 
local analysis with sub-atomic resolution of the electronic 
properties of the doped material. The STM images dis- 
play a pattern having three-fold symmetry around the 
N atoms and having a strong signal on the C atoms 
bonded to the dopant, which ab-initio simulations repro- 
duce wellJ^ STS measurements have revealed a broad 
resonant electronic state around the dopant position and 
located at an energy of 0.5 eV above the Fermi level. 8 

It is important to understand the effects of local defects 
and chemical doping on the global electronic properties 
of graphene, on its quantum transport properties and on 
its local chemical reactivity. It is the reason why the elec- 
tronic properties of graphene doped with nitrogen have 
been investigated by several groups, mainly on the basis 
of ab-initio DFT techniques &12r— The central advantage 



of the ab-initio approach is to be parameter free. A dis- 
advantage is to be restricted to periodic systems as long 
as fast Fourier transforms need to be used to link direct 
and reciprocal spaces. In most instances, doping is there- 
fore addressed in a supercell geometry that precludes the 
study of low doping concentration or disorder. In the 
case of single-wall carbon nanotubes, however, the elec- 
tronic properties of isolated defects have been studied by 
first-principle methods based on scattering theory^ For 
a nitrogen impurity in the (10,10) armchair nanotube, the 
local density of states on the N atom shows a broad peak 
centered at 0.53 eV above the Fermi level. The charge 
density associated with that quasi-bound state (donor 
level in the language of semiconductor physics) extends 
up to ~1 nm away from the defect. This means that N 
atoms located ^2 nm apart can interact, which requires 
examining with care the intrinsic validity of a supercell 
method. 

The present work, based on both ab-initio DFT and 
semi-empirical tight-binding electronic structure calcula- 
tions, aims at looking for interference effects generated by 
a distribution of N dopants in graphene as compared to 
the case of an individual impurity. We resort to different 
tight-binding parametrization strategies. The simplest 
one, based on just one adjustable parameter related to 
the defect, permits analytical calculations of the impu- 
rity levels. This model is described in Appendix A. We 
favor a more realistic approach, in which the perturba- 
tion induced by the defect is allowed to extend on the 
neighboring sites. This latter model is used to study 
the effects of disorder on the local and global electronic 
properties of doped graphene and to calculate STM im- 
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ages that we compare with available experimental data. 
The calculations and discussions are developed here for 
substitutional nitrogen. The conclusions would be qual- 
itatively the same for boron doping and for other types 
of point-like defects. 

II. METHODOLOGY 

The SIESTA package 1 ^ was used for the ab-initio DFT 
calculations. The description of the valence electrons 
was based on localized pseudo-atomic orbitals with a 
double-^ polarization (DZP). 20 Exchange-correlation ef- 
fects were handled within local density approximation 
(LDA) as proposed by Perdew and Zunger.— Core elec- 
trons were replaced by nonlocal norm-conserving pseu- 
dopotentials^ Following previous studies^^ the first 
Brillouin zone (BZ) was sampled with a 15 x 15 grid gen- 
erated according to the Monkhorst-Pack scheme^ 3 - in or- 
der to ensure a good convergence of the self-consistent 
electronic density calculations. Real-space integration 
was performed on a regular grid corresponding to a plane- 
wave cutoff around 300 Ry. All the atomic structures of 
self-supported doped graphene have been relaxed. 

As mentioned above, the DFT calculations are based 
on a code that requires a periodic system. As a conse- 
quence, a supercell scheme was adopted to handle sub- 
stitutionary doped graphene. The atomic concentration 
of dopants is then directly related to the size of the su- 
percell: 0.6 at% for a 9 x 9 supercell and 0.5 at% for a 
10 x 10 supercell. 

For the tight-binding parametrization, the following, 
well-established procedure was followed^ - — The Hamil- 
tonian included the n electrons only, the hopping pa- 
rameter between nearest-neighbor atoms was set to the 
ab-initio DFT value 70 = -2.72 eV. The same hopping 
was used between the N atom and its C first neighbors, 
which is not a severe approximation. It is indeed shown 
in Appendix A that a non-diagonal perturbation of the 
Hamiltonian can be offset by a renormalization of the N 
on-site energy, and this energy is a free parameter to be 
adjusted to the results of the ab-initio calculations. The 
C on-site energies were assumed to vary with the distance 
d to the impurity according to a Gaussian law: 

e(d) =£ C - \U\ exp(-0.5d 2 /cr 2 ) (1) 

where \U\is the depth of the potential well induced by the 
nitrogen (boron induces a potential hump, instea d 13 ! 24 ) 
and ec is the asymptotic bulk value of the on-site pa- 
rameter of carbon, ec defines the Dirac energy of pure 
graphene and has no other influence on the density of 
states (DOS) than a systematic shift of all the electronic 
levels. The standard deviation a in eq. ([1]) was set to 
0.15 nm, according to data of Ref. [24j . whereas U was 
chosen so as to best fit DFT local DOS calculated on the 
N atom in the 9x9 supercell (see below). 

In order to avoid the constraints of a periodic tight- 
binding Hamiltonian, the local DOS have been computed 



with the recursion method^ 150 pairs of recursion coef- 
ficients were calculated and extrapolated to their asymp- 
totic values related to the edges ec ± 37o of the it bands 
of pure graphene. It is worth emphasizing that this pro- 
cedure avoids to introduce any broadening of the energy 
levels. A drawback is the presence of wiggles that may 
be observed in some densities of states. They are sorts 
of Gibbs oscillations generated by the Van Hove singu- 
larities, in particular by the abrupt discontinuities of the 
graphene tt DOS at both band edges. 



III. PERIODIC DISTRIBUTION OF N 

Fig. QJa) shows DFT local densities of states on the N 
atom, on the first-neighbor atoms (CI) and on carbons 
located farther (0.5 - 0.7 nm) away (bulk), in a 9 x 9 
supercell. A broadening of 0.05 eV on a 45 x 45 grid was 
used for the calculations. The important observation is 
the existence of two peaks, located at 0.55 and 0.92 eV 
above the minimum of the density of states in the N local 
DOS4 

The present DFT calculations are in overall good 
agreement with results published for 4 x 4 , 13 ' 16 5x5^ 
7x 7p and 10 x 10 supercellsi 1 ^ The double-peak struc- 
ture in the unoccupied part of the local N DOS plotted in 
Fig. [U(a)) was not observed in previous calculations due 
to the large energy broadening that was used there,— ^ 
except in Ref. . It will be demonstrated below that the 
observed double-peak structure is a consequence of long- 
range interactions between the N impurities that were 
reproduced periodically in the graphene sheet because of 
the supercell geometry. 

Fig. [IJb) shows tight-binding densities of states of 
the same 9x9 supercell after Lorentzian broadening 
aimed at facilitating the comparison with the DFT re- 
sults. The latter are qualitatively well reproduced by 
the 7r-electron tight-binding Hamiltonian by taking \U\ 
= 4 eV in eq. ([l). The present work aims at tailor- 
ing a simple tool for exploring the effects of the de- 
fect configurations. For that reason, achieving the best 
possible fit of DFT calculations for a specific supercell 
is superfluous. By comparison, the on-site energies in 
Ref. PJ] can be approximated by a sum of two Gaussians, 
e(d) = e c - 2.95 exp(-0.5d 2 /cr 2 ) - 0.59 exp(-0.5d 2 /cr 2 ) 
(in eV), with o\ = 0.10 nm and 0-2 = 0.39 nm. The sec- 
ond term has a longer range than the single Gaussian law 
used in the present calculations, and the potential well 
at the N location is about 10% shallower than here. 

The tight-binding densities of states of the same 9x9 
superstructure, now calculated without energy broaden- 
ing, are displayed in Fig. [U(d). A substructure clearly 
emerges where the broadened local DOS only showed 
the double-peak feature discussed here above. In addi- 
tion, Fig. QJd) reveals that there are much less electronic 
states in the region between and 1.5 eV on the second- 
neighbor atoms (C2) than on the first- (CI) and third- 
neighbor (C3) carbons (see Fig [He) for the notations). 
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FIG. 1. (a) DFT and (b) tv tight-binding local densities of 
states (LDOS) of graphene doped with nitrogen at a concen- 
tration of one impurity per 9x9 supercell (162 atoms). The 
different curves represent the local DOS on the N atom, on 
the first-neighbor carbons (CI), and on a few distant carbons 
(bulk), (c) Labels assigned to the nitrogen dopant (N) and 
its first, second, third, and fourth neighbors (CI, C2, C3, and 
C4). (d) Tight-binding LDOS of the 9 x 9 superstructure 
(same as in panel b) without energy broadening, (e) DFT 
and (f) tight-binding LDOS on and around a substitutional 
N atom in graphene with 10 x 10 periodic distribution (200 
atoms per cell). The curves, shifted vertically for clarity, cor- 
respond to the atomic sites N, CI, C2, and C3 shown in panel 
(c). For all LDOS curve plotted in this figure, the zero energy 
corresponds to the Fermi level. 



The same observation as for C2 can be made for the local 
DOS on the fourth neighbors (not shown). This alterna- 
tion between lower and higher densities of states brings 
out the fact that the two sublattices of graphene are dif- 
ferently affected by a substitutional defect (see Appendix 
A). 

The reliability of the tight-binding calculations com- 
pared to DFT can be appreciated from Figs. [T|e,f) ob- 
tained for a 10 x 10 supercell (0.5 at% concentration of 
N). The good agreement between the two approaches re- 
inforces the fact that the parametrization of the tight- 
binding on-site energies used for the 9x9 supercell de- 
scribes well the doped graphene systems and does not 
require further adjustment depending on the actual con- 
centration. Electronic states localized on and near the N 
dopant are observed within 1.5 eV from the Dirac energy, 
as for the 9x9 supercell. However, the multi-peak struc- 
ture in the local N DOS of 10 x 10 (Fig.QJf)) is totally 



different from that of 9 x 9 (Fig.QJd)), despite similar N 
concentrations (0.5 at% and 0.6 at%, respectively). This 
finding is a clear indication of the importance of interfer- 
ence effects among the supercells. 

The band structures of the DFT calculations shown in 
Fig. [2] provide us with another difference between the two 
superstructures: the 10 x 10 N-doped superstructure has 
a direct band gap of ~0.05 eV at the K point of its Bril- 
louin zone, whereas the 9x9 N-doped superstructure 
has no gap. Symmetry considerations and arguments 
from perturbation theory developed in Appendix B ex- 
plain why it is so. The K and K' points of graphene 
move to the K and K' points of the folded Brillouin 
zone of a p x p supercell when p is not divisible by 3, 
whereas they are both mapped onto the T point when p 
is an integer multiple of 3~ In the latter case, the four- 
fold degeneracy of the Dirac energy at the T point of a 
supercell of pure graphene is partly lifted by the pertur- 
bation brought about by the N atoms. As demonstrated 
in Appendix B, there remain 7r and ir* bands that cross 
the Dirac energy at T. This crossing as clearly visible in 
Fig.[2ja) for the 9x9 superstructure, which therefore has 
no gap. In direct space, a doped superstructure has point 
group symmetry Dzk . The symmetry of the T point of 
the superstructure, D^h, which is the same as that of the 
K and K 1 points of pure graphene ) 27 i 28 has irreducible 
representations of dimension one and two. A twofold 
degeneracy of some electron energy bands is therefore al- 
lowed at the zone center of a doped superstructure. By 
contrast, this degeneracy is forbidden at the K and K' 
points of the folded zone because the symmetry of the- 
ses points, from D^h in pure graphene, is reduced to C^,h 
in the doped superstructure. The latter point group has 
only one-dimensional irreducible representations. As a 
consequence, the crossing of tt and tt* bands is symme- 
try forbidden at K and K' . A periodic substitution of N 
for C therefore opens a gap in these p x p supercells for 
which p is not divisible by 3, as shown in Fig.[2jb) for the 
10 x 10 superstructure. It is demonstrated in Appendix 
B that the band gap in those superstructures scales with 
p like Eg ss V/p 2 , where V is of the order of 10 eV. Other 
local gaps of the order of V/p 2 appear in the band struc- 
tures of Figs. [2]^a,b) at the edges of the folded Brillouin 
zone. They produce pseudo-gaps and sub-peaks in the 
densities of states clearly discernible in Fig. Q] for both 
the 9x9 and 10 x 10 superstructures. As a result, super- 
cells of the order of 30 x 30 at least would be necessary 
to reproduce the characteristic features of an isolated im- 
purity with an energy resolution better than 0.01 eV. 

In the tight-binding calculations, the band gap of the 
10 x 10 superstructure is located 0.14 eV below the Dirac 
point of pure graphene (ec). The Fermi level Ef of 
the doped superstructures was found to lie 0.27 eV and 
0.25 eV above Ec for the 9x9 and 10 x 10 systems, re- 
spectively, which means 0.43 eV and 0.39 eV above the 
minimum of the DOS (the corresponding DFT values are 
0.42 and 0.36 eV, respectively). The a bands below the 
Fermi energy are completely filled, the 7r bands of the 
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FIG. 2. DFT band structures of the 9 x 9 (top) and 10 x 
10 (bottom) N-doped superstructures shown along the high- 
symmetry lines of their first Brillouin zone. The zero of energy 
coincides with the Fermi level. 
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FIG. 3. Tight-binding local n densities of states on and 
around an isolated impurity in graphene. The labels cor- 
respond to the defect site (N), the first, second and third 
neighbors (CI, C2 and C3). The curves have been shifted 
vertically for clarity. The zero of energy is ec- 

IV. ISOLATED N IMPURITY 



supercell must accommodate an extra electron brought 
about by the nitrogen. The occupied local ir DOS of 
the N atom contains 1.36 electron in both the 9x9 and 
10 x 10 supercells. The remaining 0.64 excess electron is 
distributed on the surrounding carbons, of which a total 
of 0.56 sits on the three first neighbors (CI), which are 
therefore negatively charged. 

For both 9x9 and 10 x 10 doped supercells, there is 
a localized state expelled from the ir band and located 
9.35 eV below E F , to be compared with 8.65 eV in DFT 
calculations. This localized state weights more than 25% 
(0.57 state of the N local DOS, which can accommodate 
2 electrons including the spin degeneracy). The large 
value of the weight can be understood from the argu- 
ments developed in Appendix A for a simplified tight- 
binding Hamiltonian. 

We conclude here that the supercell technique would 
require very large supercells (up to 30 x 30) if it were 
desired to describe the electronic states of an isolated 
defect. In other words, the supercell technique is badly 
convergent when the size of the cell increases, particu- 
larly in low dimensions. The corresponding long range of 
a local perturbation produces interferences between the 
images of the defect generated by the periodic boundary 
conditions (see also Ref. [HI). 



We now turn our attention to the case of an isolated N 
impurity studied by tight-binding. The local DOS on and 
around an isolated N impurity are presented in Fig. [3] 
The calculations were performed on a 150x150 cell of 
graphene (45,000 atoms) with periodic boundary condi- 
tions. By restricting the number of pairs of recursion co- 
efficients to 150, the impurities located in adjacent cells 
do not feel each other, which means that the N atoms 
are virtually isolated. In the N local DOS, there is a tall 
asymmetric peak located 0.5 eV above the Dirac energy 
(ec) of graphene which coincides here with the Fermi 
level. The peak broadens and shifts to 1 eV on the first 
neighbors and moves up further to 1.5 eV on the second 
neighbors. The local DOS on the third-neighbor atoms 
reproduces the resonance peak of the impurity, with a 
smaller amplitude. The local DOS on the fourth neigh- 
bors (not shown) bears resemblance with that of pure 
graphene. It is interesting to observe in Fig. [3] that there 
remains almost nothing of the Van Hove singularities of 
graphene at ec ± 7o 011 the N and CI sites. 

The problem of an isolated impurity in graphene can be 
addressed analytically if one simplifies the Hamiltonian 
to a point-like defect by ignoring the Gaussian dereal- 
ization of the perturbation (a — > in eq. (JTJ). Details are 
presented in Appendix A. If this simple model captures 
the essential physics of the problem, it fails in produc- 
ing all the details set up by the delocalized nature of the 
perturbation. As shown in the Appendix A, perturbing 
the first-neighbor C atoms, in addition of course to the 
substitution site, leads to a better picture. 
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FIG. 4. Configurational average of the local n densities of 
states on and around substitutional N atoms in a graphene 
sheet containing 0.5 at% of nitrogen atoms randomly dis- 
tributed on the two sublattices. The labels correspond to 
the nitrogen dopant (N), the first, second and third neighbors 
(CI, C2 and C3). The curves have been shifted vertically for 
clarity. The zero of energy is ec- 
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FIG. 5. Configurational average of the local it densities of 
states on and around substitutional N atoms in a graphene 
sheet containing 0.5 at% of nitrogen atoms randomly dis- 
tributed on the same sublattice. The labels correspond to 
the nitrogen dopant (N), the first, second and third neighbors 
(CI, C2 and C3). The curves have been shifted vertically for 
clarity. The zero of energy is ec- 



V. RANDOMLY DISTRIBUTED N ATOMS 

Experimentally, nitrogen doping of graphene leads to 
randomly distributed substitutional sites*^ with per- 
haps some preference of the dopant atoms to sit on the 
same sublattice, at least locally^ The calculations illus- 
trated in the present Section were performed for nitro- 
gen distributed randomly with an atomic concentration 
of 0.5%, identical to that realized with the 10 x 10 su- 
perstructure. The selection of the substitutional sites 
was constrained by the requirement that the distance be- 
tween two nitrogens remains larger than 12a = 1.8 nm 
(see eq. |T])). The reason for that was to avoid any over- 
lapping of the potential wells generated by the dopants. 

Fig. @] shows a configurational averaging of the local 
densities of states when the N atoms have equal proba- 
bilities to sit on one or the other of the graphene sublat- 
tices. This case will be referred to as "unpolarized" (or 
"uncompensated")^ By configurational averaging, it is 
meant the arithmetic average of local DOS on 50 nitro- 
gens selected randomly — among the 225 dopant atoms 
that the 150x150 sample box contains — and the arith- 
metic average of local DOS on 150 randomly-selected CI, 
C2 and C3 sites. The N, CI, C2, and C3 average local 
DOS for the unpolarized disordered distribution all have 
a remarkable similarity with the ones obtained for an 
isolated impurity (Fig. [3]) . The result are in agreement 
with N local DOS for randomly distributed impurities 
obtained by Lherbier et ali^ with a small broadening. 



When the N atoms are put, still randomly, on the same 
sublattice ("polarized" case"), the multi-peak structure 
characteristics of the superlattice (see Fig. [TJe,f)) reap- 
pear on the average local DOS shown in Fig. [5j Interest- 
ingly, the local DOS on the second-neighbor atoms (C2), 
which belongs to the same sublattice as the N atom they 
surround, are virtually the same in Figs. U and [SJ in- 
sensitive to the polarization of the N distribution. It is 
also interesting to remark that, like with the 10 x 10 su- 
perstructure, there is a tiny gap of states 0.1 eV below 
the Dirac energy. The gap looks narrower in case of the 
unpolarized distribution compared to the polarized one. 
In the polarized case, the appearance of a gap is natural 
since the average diagonal potential breaks the symme- 
try between the two sublattices. Similar effects have been 
observed in the case of vacancies^ 

Averaging the local densities of states makes visible 
the distinction between unpolarized and polarized distri- 
butions of N. However, for a given configuration of the 
dopant atoms, the local DOS varies from site to site. 
How much that variation can be is illustrated in Fig. [51 
which compares local DOS calculated on three different 
N atoms for both the unpolarized and the polarized dis- 
tributions. The variations from site to site are so impor- 
tant that the local DOS is not a reliable indicator of the 
global distribution of the nitrogen atoms among the two 
sublattices. If it is true that the N local DOS displayed 
in Fig. [B] are less peaked that the N local DOS of the 
10 x 10 superstructure (Fig [Ue,f)), identifying which is 
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FIG. 6. Local DOS on three N atoms randomly selected in 
graphene doped at concentration 0.5 at%. The unpolarized 
case (left) refers to dopants located on the two sublattices, 
whereas the dopants all sit on the same sublattice in the po- 
larized case (right). The zero of energy is ec- 



which would be difficult on the experimental side if one 
had only STS spectra to deal with. 

Finally, it may be important to remark that the min- 
imum of the 7r density of states around the N dopants 
does not coincide with the atomic level ec of carbon in 
graphene (the Dirac point energy; see also Refs. [3(3, HH, 
l32j|). Depending on the concentration, the DOS mini- 
mum shifts slightly below ec (see e.g. Figs. [4] and [5]). 



VI. STM IMAGES 

A reliable information on the electronic structure of 
doped graphene can be obtained by STM imaging. A 
tip polarized negatively compared to the graphene layer 
probes the unoccupied states where a N substitutional 
impurity and the adjacent carbon atoms have peaks in 
their local DOS. The increase of electronic density of 
states, compared to graphene, in this energy region gives 
rise to a bright triangular spot in the STM image , 5 i 8 ] 12 i 17 
with two possible orientations with respect to the honey- 
comb lattice depending on the sublattice on which the N 
dopant sits. These two orientations are rotated by 180° 
from each other, as observed experimentally.—^ 

Fig. [7] is a tight-binding STM imaged computed for 
graphene with a single N impurity. The prominent tri- 
angular arrow-head at the center of the image is located 
on the N, CI and C3 sites. The second-neighbor carbons 
(C2) atoms do not participate much to the STM signa- 
ture of the dopant. This is so because the local DOS 
in the energy window probed by the STM current (0 to 
0.5 eV) is smaller on the C2 atoms than on the N, CI 
and C3 atoms (see Fig. |3J). 
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FIG. 7. Tight-binding STM constant-current image of 
graphene with an isolated N impurity located near the center, 
computed for a negatively-polarized tip (Vtip = -0.5 V). The 
vertical scale is the tip height (A) above sample. The image 
size is 2 nm along both sides. 



Fig. [5] shows two STM images computed for graphene 
doped with N at 0.5 at% concentration, the dopants be- 
ing randomly distributed on the two sublattices. In the 
configuration shown in image (a), there are two impuri- 
ties located 1.9 nm apart sitting on the same sublattice. 
The two nitrogens captured in images (b) do not sit on 
the same sublattice, which explains why the related tri- 
angular STM patterns are oriented differently. The scale 
used for the STM signal (tip height at constant current) 
is identical for both images. For improving the contrast, 
the scale has been saturated below 4.0 A and above 5.8 A. 
Two dopants interfere more when they are located on the 
two sublattices than when they sit on the same sublat- 
tice. 

In tight-binding STM theory^ there is an empirical 
parameter coupling the tip apex to the sample atoms, 
which was chosen independent of the chemical nature of 
the probed atoms. In other words, the STM calculations 
differentiate a N atom from the C atoms only through in- 
trinsic effects that the former has on the electronic struc- 
ture of the doped graphene. In DFT calculations, the 
standard way to generate an STM image is via Tersoff- 
Hamann's theory^ In addition to local variations of the 
DOS, the tunnel current depends on the spatial exten- 
sion of the 7r orbitals in the direction perpendicular to 
the atomic plane. For doped graphene, the 2p z orbital of 
N decreases more rapidly than the C 2p z orbital does^ 
which means that the contrast of the computed image 
depends on the tip-sample distance. This dependence 
is illustrated in Figs. |H{a) and (b) that represent DFT 
constant-current STM images computed for a small cur- 
rent (large tip-sample distance) and for a large current 
(small tip-sample distance), respectively. The local DOS 
is larger on top of the N atom than on the top of all 
the other atoms, but it decreases much more rapidly in 
the normal direction. This explains the small brightness 
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FIG. 8. Constant-current STM images computed for two con- 
figurations of graphene doped with 0.5 at% N. The distribu- 
tion of the N atoms is random. The two nitrogen atoms visible 
in image (a) lie on the same sublattice, which is not the case 
in image (b). In both cases, the STM tip is polarized at -0.5 V 
with respect to the sample. Each image is a square of 4 nm 
edge. 



of the N site for the large tip-sample distance. Experi- 
mentally, variations of the STM contrast have also been 
observed as a function of the current /voltage conditions, 8 
where image of the dopant atom protrudes in some occa- 
sions and does not in other occasions. It is tempting to 
attribute this observation to variations of the tip-sample 
distances depending on the current setpoint, as in Fig. ® 
However, the experimental distance is much larger than 
the one that can be achieved numerically in DFT, due 
to the fast decay of the localized basis set when moving 
away from the sample surface. 



VII. CONCLUSIONS 

Tight-binding and DFT electronic calculations reveal 
that the perturbation induced by N dopants in graphene 
is long ranged. This is not a surprise when one re- 
members that the Green function of the free-electron 
problem in two dimensions is the Bessel-Hankel func- 
tion of zeroth order (i/4).Hq (As|f* — f*\), which slowly 
decreases like 1/ykd with the distance d = \r — r*], where 




FIG. 9. DFT calculation of constant-current STM images of 
10 x 10 N-doped graphene for two levels of the tunnel current: 
(a) high and (b) low (see text). The tip bias polarization is 
-0.5 V, the vertical scale is the tip height (A) above sample. 
Each image represents a square of 4 nm edge centered on a N 
dopant. 



k = \J h 2 /2mE. For the 7r-electron tight-binding Hamil- 
tonian on the honeycomb lattice, the Green function el- 
ements G° m can also be approximated, at low energy, 
by Hankel functions of Sk\f m — f n \, now of order zero or 
one, depending on whether they refer to sites m and n lo- 
cated on the same sublattice or not<^ In this regime, the 
energy E is linearly related to 8k by the dispersion rela- 
tion E = hvpSk, where vf is the Fermi velocity, and the 
Hankel functions are also multiplied by the energy and 
by periodic functions describing -\/3 x -\/3 modulations 
related to the Dirac wave vectors K and K'. For large 
separation distances, the Green function elements decay 
slowly, like 1/ ' \JE \f m — r n \/hvF, except for very small 
excitation energies where the Hankel functions diverge. 
A more detailed discussion will be given elsewhere^ 

The scattering formalism developed in Appendix A for 
the case of a point-like defect emphasizes the central role 
of the Green function in the understanding of the per- 
turbation induced by the defect. Hence the long-range 
interaction between dopants in graphene. It is impor- 
tant to realize that this long-range interaction is not the 



8 



consequence of using an extended perturbation of the on- 
site energies around the defect (eq. JT])). Things actually 
go the other way round: a defect perturbs the crystal 
potential far away, because of the slow decay of the scat- 
tering it produces. As a consequence, the parameters of 
the tight-binding Hamiltonian that mimic ab-initio cal- 
culations are affected in a sizable neighborhood of the 
defect^ 

This long-range interaction has several effects. First, 
it produces interferences between the duplicates of the 
defect generated by periodic boundary conditions when 
a supercell approach is used. The latter must therefore 
be used with caution. Second, substitutional impurities 
dispersed in graphene cannot be treated as independent 
defects as soon as the distance between two of them be- 
comes too small, about 2 nm in the case of nitrogen. In- 
creasing the strength of the local potential also increases 
the range of the perturbation it produces. This unusual 
behavior is actually related to the properties of the Green 
functions at low energy (see Appendix A and Refs. [HI, 
US Hl|). As mentioned in Appendix A, the deepest de- 
fect is the vacancy. Nitrogen, with its 10 eV effective 
perturbation parameter, is not a shallow defect. Boron 
would not be a shallow impurity either. B substitution 
can be described by a potential hump, instead of a well, 
with a similar \U\ as nitrogen and a similar long-range 
perturbationj 13 i 24 In the case of more complex defects, 
such as N plus vacancy^ P plus N^S or it would 

be more difficult to define local parameters and to assess 
the importance of long-range interactions along the same 
way as in this paper. 

As demonstrated experimentally, the partition of the 
honeycomb network in two sublattices has subtle effects 
on STM image of graphene with substitutional impuri- 
ties, as exemplified by Fig.[5J The STM image of a dopant 
in graphene may be influenced by the proximity of an- 
other dopant. This is a direct consequence of the defect 
interactions mediated by graphene. What happens when 
two impurities come very close to each other remains to 
be clarified. A new parametrization of the tight-binding 
Hamiltonian would be required if one had to obtain, 
from calculations, the STM topography of neighboring 
N dopants. 
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APPENDIX A. THE IMPURITY PROBLEM IN 
GRAPHENE 



In this Appendix, a simple but illustrative perturba- 
tion model is developed. It is assumed that the on-site 
energies all take the unperturbed bulk value Ec = 0, 
except on the nitrogen atom where Sn — U. The tight- 



binding Hamiltonian can be set in the form: 

H = H Q + V = -Y J Wnm{m\ + \{))U{{)\ . (2) 

n,m 

The states \n) denote tt orbitals centered on sites n: 
(^(r — n) = (r\n). H is the usual tight-binding Hamil- 
tonian where only hopping integrals t nm = 70 between 
first neighbors n and m are kept, and V is the localized 
potential of the N atom at site 0. The local density of 
states n(r, E) is given by: 

n(r,E) = y]rinm(E)(pn(r - n)ip^(r - m) , (3) 

n, m 

where n nm (E) is obtained from the Green function or 
resolvent G(z) = (z — 



n nm {E) 



— lim Im (n\G(E + ie)\m) . (4) 

7T £->0 



G(z) can then be calculated in terms of the unperturbed 
Green function G°(z) = (z — H )^ 1 ; this is the so-called 
Koster-Slater-Lifshitz problem. If the matrix elements 
(n\G(E + ie)\m) are written G nm {z), we have: 



G, 



L, n0 TU 0m ' 



u 



1 - UG° 00 



(5) 



In particular, on the impurity site, we have Goo = 
Goo/i 1 ~ UG% )j%&. and the local density of states on 
the N site, um{E) = noo(E) is given by: 



n N {E) 



n°{E) 



(1 - UF°{E)) 2 +Tr 2 U 2 n°(EY 



(6) 



where F°(E) is the real part of the Green function 
Gq (E + ie), i.e. the Hilbert transform of the unper- 
turbed density of states of graphene n°(E). 

Close to the Dirac point, n°(E) ~ \E\ and F°(E) ~ 
£71n|.E|, and provided U is large enough, un{E) shows 
a resonant behavior around the energy E r such that 
(l-UF°(E r )) = (see the open-circle symbol in Fig. [T0| . 
This has been discussed in many places i 30 i 32 ' 35 i 37 ' 44 ~ — It 
turns out that a value for U of the order of —10 eV re- 
produces a resonance at the position obtained with the 
previous model (compare Fig. [TlTa) with Fig. [3]) instead 
of the value U = —4 eV in the previous approach. How- 
ever, the shapes of the local DOS differ to some extent be- 
tween the local-perturbation model developed in this Ap- 
pendix and the more delocalized potential well considered 
throughout this paper. The resonance here (Fig. [Ufa)) 
is weaker on the N site than on the sites CI and C3 lo- 
cated on the other sublattice, which can be explained 
from the symmetry properties of the unperturbed Green 
functions: G nm {z) is and odd or even function of z de- 
pending on whether the sites n and m belonging to the 
same sublattice of the graphene structure or not. As a 
consequence, for example, Gq Q (z = 0) = whereas the 
real part of Gq u (z = 0) does not vanish when n is a first 
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FIG. 10. Real (F°(E)) and imaginary {n°(E)) parts of the 
Green function of the graphene n tight-binding Hamiltonian 
along the real axis, in units of |-yo | - The intersection of F°(E) 
with the reciprocal of the perturbed level 1/U (dashed line) in 
the vicinity of the Dirac point gives the energy E r at which 
the local DOS of the impurity has a resonance (open circle 
symbol). There is another intersection (filled symbol) giving 
rise to a localized state below the 7r-electron density of states. 



neighbor, which has a big impact on the local densities 
of states (see eq. ((S|))i^ Also the Van Hove singularities 
at 2.7 eV are still present here. 

In addition to the resonant state located above the 
Dirac point, Fig. [TUJ reveals the existence of a localized 
state lying below the 7r-electron DOS, where the on- 
site perturbation level 1/U intersects the curve F°(E) at 
— 4.5|7o| (solid square symbol). The weight (residue) of 
this localized state is inversely proportional to the slope 
of F°(E) at the intersection. This is the localized state 
found at 9.35 eV below Ep in tight-binding and 8.65 eV 
below Ep in DFT calculations for the doped superstruc- 
tures (Section UlI)) . 

This one-parameter model can be generalized to other 
type of local defects or chemical doping. B doping, with 
a positive U, will yield results symmetrical to those ob- 
tained for N: it can indeed be anticipated from Fig. QJJ 
that the resonance state will now appear near the top 
of the occupied states, as confirmed by DFT calcula- 
tions The limit U — > oo corresponds to the intro- 
duction of a vacancy instead of a substitutional impurity. 
Because of the vanishing of the density of states and of 
the logarithmic divergence of F°(E), this limit is singular 
and must been studied with special carei 30 i 32 ' 45 ~ — The 
resonant state becomes very narrow, but does not become 
a genuine bound state. The total density of states varies 



as l/([ln |_E|] 2 |_E|) and the perturbation of the electronic 
density is concentrated on the first neighbors and on the 
sites belonging to the sublattice different from that of the 
vacancy^ 

The simple model developed here above allows one to 
address the question of the role of N-C hopping inter- 
actions on the local DOS on the impurity. When the 
N-C hopping takes a value 7 different from the one 70 
between the C atoms, the Green function element on the 
perturbed site can be set in the form 



Goo( z ) 



1 



z-[/-(7 2 /7o 2 )£°M 



(7) 



where S°(z) is the self energy of the unperturbed 
graphene. In the notations of the above formalism, the 
self energy can be identified to £°(z) = z — 1/Gq (z). 
Inserting this expression in the right-hand side of eq. ([7]) 
leads to a generalization of eq. ([6]) valid for 7 7^ 7o- In 
particular, the resonance condition in the vicinity of E 
= becomes 



7o^+(7 2 -7o 2 )£. 



F°(E r ) 



(8) 



By comparison with the situation depicted in Fig. 1101 
the intersection of F°(E) needs now to be search for 
with a curve whose ordinate and slope at the origin are 
-l 2 l(ll\U\) and 7 2 (7o 2 " 7 2 )/(7 2 ^) 2 , respectively. If 
I7I < 1 7o I ! which should be the case here since the N atom 
is smaller than the C one, the ordinate at the origin of the 
curve moves up and its slopes becomes positive. These 
two effects pull the resonant energy E r closer to the origin 
compared with the simplest situation where I7I = |7o|. 
However, if the hopping perturbation is small, renormal- 
izing the on-site level U to a larger absolute value U e g 
would produce the same effect. This conclusion validates 
the approach used so far to not modifying the N-C hop- 
ping. 

We have finally considered an intermediate model be- 
tween the Gaussian distribution of the on-site energies 
and the one-parameter model just described. We took 
Uq on the N site and a second perturbation potential 
Ui on the first neighbors, with values fixed by eq. ([lj, 
U a = -4 eV; Ui = -2.57 eV. The agreement with the 
full Gaussian model is much better than with the one 
parameter model (Fig. [TlTb)). It is fairly remarkable to 
realize that the intensity of the resonance state on the 
N site increases when the potential is delocalized on the 
first neighbors. 



APPENDIX B. BAND GAP OPENING IN 
DOPED SUPERSTRUCTURES 

The following perturbation theory is based on the same 
model Hamiltonian as in eq. ([2]) , except that the potential 
V refers now to the perturbation brought about by the 
periodic array of nitrogen atoms in substitution for C in 
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0.25 




FIG. 11. Local densities of states deduced from the sim- 
ple tight-binding model for an isolated nitrogen impurity in 
graphene. (a) On-site perturbation energy U = -10 eV on the 
N site only; (b) on-site perturbation on the N (-4 eV) and the 
three CI sites (-2.57 eV). 



a p x p supercell of graphene. These atoms occupy one 
among the 2N = 2p 2 sublattices of the supercell, here 
denoted sublattice 1: 



H = H° + V = - ^2\n)t, 



(9) 



n£l 



When V — 0, there are four states with zero energy 
(Dirac energy of the unperturbed graphene): \K % ) and 
\K n ), where i = A,B denote the sublattices of the 
graphene structure and \K l ), \K H ) are the corresponding 
Bloch functions: 



\ K A(B)\ 



1 



E 

ieA(B) 



iK. 



(10) 



where K.n in the exponential represents the dot product 
of the wave vector of the K point of graphene and the 
position vector of the site n in real space. 



It will be assumed that the sublattice 1 is an A sub- 
lattice. Then, the states \K B ) and \K' B ), having no am- 
plitude on sublattice 1, are eigenstates of zero energy for 
any value of the on-site perturbation U. Very close to 
zero energy, furthermore, only the ((\K A ), \K' A ))) sub- 
space needs to be considered. In particular, 

(K A \V\K A ) = (K' A \V\K' A ) = I ]T 1 = U/p 2 (11) 



nel 



(K A \V\K' A ) = i £exp[i(** - K).n] = U /p 2 S K > -liW) 

nel 

where the Kronecker symbol imposes that K' — K be a 
vector of the reciprocal lattice of the supercell, i.e. it is 
one when p = 3g is a multiple integer of 3 and is if p is 
not divisible by 3. 

As a consequence, in the case p = 3g, the four states 
with zero energy are mapped onto the zone center. Two 
of them, with non-zero amplitudes on the B sites, remain 
degenerate with zero energy; the other two states have 
energy E = and E = —2\U\/p 2 to lowest order in 
perturbation. There are then two plus one states at zero 
energy, but this is an accidental degeneracy due to the 
simplified tight-binding model. When p is not divisible by 
3, the degeneracy at the K and K' points of the folded 
zone is lifted: at each point, one state remains at zero 
energy (\K B ) or \K' B ), respectively) whereas the other 
state moves down at energy — \U\/p 2 to lowest order in 
perturbation. 

All these features are clearly visible in the band struc- 
tures of the 9x9 and 10x10 superlattices shown in Fig. [2] 
Since a realistic value for \U\ is about 10 eV, the gap for 
the 10 x 10 case should be of the order of 10/p 2 = 0.1 eV, 
not far from the actual value (see Section III). For the 
9x9 superstructure, the expected crossing between the 
two linear ir and ir* bands is observed at the T point close 
to two parabolic branches separated from each other by 
a gap about twice this value. 
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